In [1]:
%pylab inline
import matplotlib.pyplot as plt

from IPython.core.display import HTML
Populating the interactive namespace from numpy and matplotlib

For details see Chapter 22. in Jenő Sólyom: Fundamentals of the Physics of Solids, Volume 2: Electronic Properties (A modern szilárdtes-tfizika alpajai II. Fémek, félvezetők, szupravezetők).

In [2]:
# Az abra kimentesehez az alabbiakat a plt.show()  ele kell tenni!!! 

#savefig('fig_rainbow_p2_1ray.pdf');  # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps');  # Abra kimentese

# Abra es fontmeretek
xfig_meret= 9   #    12 nagy abrahoz
yfig_meret= 6    #   12 nagy abrahoz
xyticks_meret= 15  #  20 nagy abrahoz
xylabel_meret= 30  #  30 nagy abrahoz
legend_meret= 21   #  30 nagy abrahoz
In [3]:
def nmax_func(x):
    nn=floor(x)
    if(x-nn <= 1/2):
        nmax=nn-1
    else:
        nmax=nn
    
    #print ("nmax+1/2 < x < nmax+3/2", nmax+1/2, x, nmax+3/2)
    return (nmax)
In [4]:
def mu_func(x):
    
    nmax = nmax_func(x)
    nlist=arange(nmax+1)
     
    tmp=sum(sqrt(x-(i+1/2)) for i in nlist)

    #print("sum=",res)
    #print(nmax,nn,xx)
    
    res=(3/2*tmp)**(2/3)
    
    return res
In [5]:
def dos(x):
    nmax = nmax_func(x)
    nlist=arange(nmax+1)
        
    res=0.5*sum(1/sqrt(x-(i+1/2)) for i in nlist)
    return res
In [6]:
def dos_Fermi(x):
    # x ---> hbar omega_c/E_F 
    
    nmax = floor(1/x-1/2)
    nlist=arange(nmax+1)

    #print("nmax= ",nmax)
    res=0.5*sqrt(x)*sum(1/sqrt(1/x-(i+1/2)) for i in nlist)
    
    return res
In [7]:
#  energy in 2D
def en_2D(B):
    n=floor(1/B)
    res=2*B*(n+1/2-1/2*n*(n+1)*B)
    
    return res
In [8]:
#  M in 2D
def M_2D(B):
    n=floor(1/B)
    res=2*B*(n+1/2-1/2*n*(n+1)*B)
    
    return res
In [9]:
#  energy in 3D
def En_Landau_dia(x):
    # energy in units of N_e * E_F
    
    # x ---> hbar omega_c/E_F
    # x <= 2
    
    nmax = floor(1/x-1/2)
    nlist=arange(nmax+1)

    #print("nmax= ",nmax)
    res=1-x/2*sum((1-(i+1/2)*x)**(3/2) for i in nlist)
    
    return res
In [10]:
#  M in 3D
def M_Landau_dia(x):
    # x ---> hbar omega_c/E_F 
    
    nmax = floor(1/x-1/2)
    nlist=arange(nmax+1)

    #print("nmax= ",nmax)
    res=1/2*sum((1-(i+1/2)*x)**(3/2) for i in nlist)
    
    return res
In [11]:
# chi in 3D, susceptibility
def chi_Landau_dia(x):
    # x ---> hbar omega_c/E_F 
    
    nmax = floor(1/x-1/2)
    nlist=arange(nmax+1)

    #print("nmax= ",nmax)
    res=-3/2*sum((i+1/2)*(1-(i+1/2)*x)**(1/2) for i in nlist)
    
    return res
In [12]:
def M_osc(B,T,params):
    
    # Solyom konyv II. (22.4.86)
    # B---> hbar omega_c/E_F
    # T---> k_B T/hbar omega_c
    # mcs ---> m^*/m)e
    # ge ---> g_e
    ge,mcs,lmax = params
    llist=arange(1,lmax+1)

    sum = 0;
    for l in llist:
        x=1/sqrt(l)*(-1)**(l+1)*cos(pi*l*ge/2*mcs)*sin(pi/4-2*pi*l/B)/(sinh(2*pi*pi*l*T))**2    
        sum = sum + x
    
    res=(ge/2)**2-1/3*(1/mcs)**2 + 2*T*sqrt(2/B)*sum
    
    return res

Density of states (DOS) $\varrho(E)$ in 3D

$\frac{\varrho(E)}{\varrho_0(E)} = \frac{1}{2}\, \sqrt{\frac{\hbar \omega_c}{E}}

\sum{n=0}^{n{\mathrm{max}}}\, \frac{1}{\sqrt{\frac{E}{\hbar \omega_c}-(n+\frac{1}{2})}}$,

where $\omega_c = \frac{eB}{m}$, $(n_{\mathrm{max}}+1/2)\,\hbar \omega_c \le E$ and $\varrho_0(E)$ is the DOS with zero magnetic field.

In [13]:
Npoints=100;
en_max=7;
enlist = linspace(0,en_max,Npoints)
gyoken = sqrt(enlist)
 
dd = [];

for i in enlist:
    dd.append(dos(i))

plot(enlist,dd,label=r'$\varrho(E)$',color='r')  
plot(enlist,gyoken,color='b')  

#legend(loc='upper left',fontsize=legend_meret)

title("DOS",fontsize=20)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\frac{E}{\hbar \omega_c}$',fontsize=xylabel_meret)
ylabel(r'$\varrho(E)$',fontsize=xylabel_meret)

xlim(0,en_max*1.)
ylim(0,en_max*1.)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese

B dependence of the DOS at the Fermi energy

$\frac{\varrho(E_\mathrm{F})}{\varrho0(E\mathrm{F})} = \frac{1}{2}\, \sqrt{\frac{\hbar \omegac}{E\mathrm{F}}}

\sum{n=0}^{n{\mathrm{max}}}\, \frac{1}{\sqrt{\frac{E_\mathrm{F}}{\hbar \omega_c}-(n+\frac{1}{2})}}$, where $\omega_c = \frac{eB}{m}$ and $\varrho0(E\mathrm{F})$ is the DOS at the Fermi energy with zero magnetic field.

In [14]:
Npoints=200;
B_max=1;
Blist = linspace(0.1,B_max,Npoints)
 
dd = [];

for i in Blist:
    dd.append(dos_Fermi(i))

plot(Blist,dd,label=r'$\varrho(E_\mathrm{F})$',color='r')  
#plot(enlist,gyoken,color='b')  

#legend(loc='upper left',fontsize=legend_meret)

title("DOS at the Fermi energy",fontsize=20)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\frac{\hbar \omega_c}{E_\mathrm{F}} \sim B$',fontsize=xylabel_meret)
ylabel(r'$\varrho(E_\mathrm{F})$',fontsize=xylabel_meret)

xlim(0,B_max*1.)
#ylim(0,1.5)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese

1/B dependence of the DOS at the Fermi energy

equidistant peaks

In [15]:
Npoints=300;
B_max=5;
Blist = linspace(0.01,B_max,Npoints)
 
dd = [];

for i in Blist:
    dd.append(dos_Fermi(1/i))

plot(Blist,dd,label=r'$\varrho(E_\mathrm{F})$',color='r')  
#plot(enlist,gyoken,color='b')  

#legend(loc='upper left',fontsize=legend_meret)

title("DOS at the Fermi energy",fontsize=20)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\frac{E_\mathrm{F}}{\hbar \omega_c} \sim 1/B$',fontsize=xylabel_meret)
ylabel(r'$\varrho(E_\mathrm{F})$',fontsize=xylabel_meret)

xlim(0,B_max*1.)
ylim(0,5)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese

B dependence of the chemical potential $\mu(B):$

$\eta = \frac{\mu(B)}{\hbar \omega_c}$ and $\eta_0 = \frac{E_\mathrm{F}}{\hbar \omega_c}$

$\sum_{n=0}^{n_{\mathrm{max}}}\, \sqrt{\eta-(n+\frac{1}{2})} = \frac{2}{3}\, \eta_0^{3/2}$

where $n_{\mathrm{max}} +1/2 < \eta < n_{\mathrm{max}} +3/2$

For each values of $\eta$ we calculate $\eta_0$ from the above equation then $\eta$ is plotted as a function of $\eta_0$.

In [16]:
ss=[1.5,2.5,3.5]

kk=[];
ll=[];
pp=[];
dd=[];

for i in ss:
    ll.append(nmax_func(i))
    kk.append(mu_func(i))
    pp.append(i/mu_func(i))
    dd.append(dos(i))

print(ll)
print (kk)
print(pp)
print(dd)
[0.0, 1.0, 2.0]
[1.3103706971044482, 2.3581854906485358, 3.3819368911950205]
[1.1447142425533319, 1.0601371308210632, 1.0349099089082237]
[0.5, 0.85355339059327373, 1.1422285251880866]
In [17]:
Npoints=100;
eta_max=5;
eta = linspace(0,eta_max,Npoints)
 
eta_0 = [];

for i in eta:
    eta_0.append(mu_func(i))

plot(eta_0,eta,label=r'$\eta = \frac{\mu(B)}{\hbar \omega_c}$',color='r')  
plot(eta,eta,color='b')  

legend(loc='upper left',fontsize=legend_meret)

title(r'$\mu(B)$',fontsize=xylabel_meret*0.8)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\eta_0 = \frac{E_\mathrm{F}}{\hbar \omega_c}$',fontsize=xylabel_meret)
ylabel(r'$\eta$',fontsize=xylabel_meret)

xlim(0,eta_max*1.1)
ylim(0,eta_max*1.1)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese

B dependence of the chemical potential $\frac{\mu(B)}{E_\mathrm{F}}$

$\frac{\mu(B)}{E_\mathrm{F}}= \frac{\eta}{\eta_0}$

The result shows that $\mu(B)$ is weakly depend on B and practically equals to $E_\mathrm{F}$ for not too large magnetic field.

In [18]:
Npoints=200;
eta_max=5;
eta = linspace(0.01,eta_max,Npoints)
unit=linspace(1,1,Npoints)

etap = [];
eta_0 = [];

for i in eta:
    eta_0.append(mu_func(i))
    etap.append(i/mu_func(i))

plot(eta_0,etap,label=r'$\frac{\mu(B)}{E_\mathrm{F}}$',color='r') 
plot(eta_0,unit,color='b')

legend(loc='upper right',fontsize=legend_meret)

title(r'$\mu(B)$',fontsize=xylabel_meret)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\eta_0 = \frac{E_\mathrm{F}}{\hbar \omega_c}$',fontsize=xylabel_meret)
ylabel(r'$\frac{\mu(B)}{E_\mathrm{F}}$',fontsize=xylabel_meret)

xlim(0.,eta_max)
ylim(0,2)
grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese
/usr/local/lib/python3.4/dist-packages/ipykernel/__main__.py:11: RuntimeWarning: divide by zero encountered in double_scalars

Energy in 2D as a function of B

In [19]:
Npoints=300;
B_max=1.1;
Blist = linspace(0,B_max,Npoints)

dd = [];

for i in Blist:
    dd.append(en_2D(i))

plot(Blist,dd,label=r'$\varrho(\epsilon)$',color='r')  
#plot(enlist,gyoken,color='b')  

#legend(loc='upper left',fontsize=legend_meret)

title(r'$E_0$',fontsize=20)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\frac{B}{B_0}$',fontsize=xylabel_meret)
ylabel(r'$E_0$',fontsize=xylabel_meret)

xlim(0,1.1)
#ylim(0,eta_max*1.)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese
/usr/local/lib/python3.4/dist-packages/ipykernel/__main__.py:3: RuntimeWarning: divide by zero encountered in double_scalars
  app.launch_new_instance()
/usr/local/lib/python3.4/dist-packages/ipykernel/__main__.py:4: RuntimeWarning: invalid value encountered in double_scalars

Energy in 2D as a function 1/B

In [20]:
Npoints=400;
B_max=6.27;
Blist = linspace(0.5,B_max,Npoints)

dd = [];

for i in Blist:
    dd.append(en_2D(1/i))

plot(Blist,dd,label=r'$\varrho(\epsilon)$',color='r')  
#plot(enlist,gyoken,color='b')  

#legend(loc='upper left',fontsize=legend_meret)

title(r'$E_0$',fontsize=20)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\frac{B_0}{B}\sim 1/B$',fontsize=xylabel_meret)
ylabel(r'$E_0$',fontsize=xylabel_meret)

xlim(0,B_max)
ylim(1,1.15)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese

Total energy in 3D as a function of the magnetic field at $T=0$

In [21]:
En_Landau_dia(0.27)
Out[21]:
0.80267936944445906
In [22]:
Npoints=200;
B_max=2.;  # hbar omega_c/E_F -->  B_max <= 2
Blist = linspace(0.01,B_max,Npoints)
 
dd = [];

for i in Blist:
    dd.append(En_Landau_dia(i))

plot(Blist,dd,label=r'$E_0(B)$',color='r')  

#plot(enlist,gyoken,color='b')  

#legend(loc='upper left',fontsize=legend_meret)

title("Total energy $E_0(B)$ at $T=0$",fontsize=20)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\frac{\hbar \omega_c}{E_\mathrm{F}} \sim B$',fontsize=xylabel_meret)
ylabel(r'$E_0(B)$',fontsize=xylabel_meret)

xlim(0,B_max*1.)
#ylim(0,0.5)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese

Landau diamagnetization in 3D as a function of $B$ at $T=0$

In [23]:
M_Landau_dia(0.27)
Out[23]:
0.73081715020570726
In [24]:
Npoints=200;
B_max=1;  # hbar omega_c/E_F -->  B_max <= 2
Blist = linspace(0.1,B_max,Npoints)
 
ddm = [];

for i in Blist:
    ddm.append(M_Landau_dia(i))

plot(Blist,ddm,label=r'$M_(B)$',color='b')  

#plot(enlist,gyoken,color='b')  

#legend(loc='upper left',fontsize=legend_meret)

title("Total magnetization $M(B)$ at $T=0$",fontsize=20)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\frac{\hbar \omega_c}{E_\mathrm{F}} \sim B$',fontsize=xylabel_meret)
ylabel(r'$M(B)$',fontsize=xylabel_meret)

xlim(0,B_max*1.)
#\ylim(0,7)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese

Landau diamagnetic contribution to the susceptibility $\chi(B)$ at $T=0$ and in 3D

In [25]:
Npoints=200;
B_max=1;  # hbar omega_c/E_F -->  B_max <= 2
Blist = linspace(0.1,B_max,Npoints)
 
ddm = [];

for i in Blist:
    ddm.append(chi_Landau_dia(i))

plot(Blist,ddm,label=r'$M_(B)$',color='b')  

#plot(enlist,gyoken,color='b')  

#legend(loc='upper left',fontsize=legend_meret)

title("Landau dia, susceptibility $\chi(B)$ at $T=0$",fontsize=20)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\frac{\hbar \omega_c}{E_\mathrm{F}} \sim B$',fontsize=xylabel_meret)
ylabel(r'$\chi(B)$',fontsize=xylabel_meret)

xlim(0,B_max*1.)
#\ylim(0,7)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese

Magnetization oscillation with B at finite temperature T and in 3D

In [26]:
B=0.17;
T=0.027;

ge=2;
mcs=2;
lmax =10;

print("B, T, ge, mcs, lmax = ",B, T, ge, mcs, lmax)

params=[ge, mcs, lmax]

print(params)
M_osc(B,T,params)
B, T, ge, mcs, lmax =  0.17 0.027 2 2 10
[2, 2, 10]
Out[26]:
1.4526382561919893
In [27]:
T=0.01;
ge=2;
mcs=1;
lmax =10;

print("T, ge, mcs, lmax = ", T, ge, mcs, lmax)

params=[ge, mcs, lmax]

Npoints=200;
B_max=0.5;  # hbar omega_c/E_F -->  B_max <= 2
Blist = linspace(0.05,B_max,Npoints)
 
ddm = [];

for i in Blist:
    ddm.append(M_osc(i,T,params))

plot(Blist,ddm,label=r'$M_(B)$',color='b')  

#plot(enlist,gyoken,color='b')  

#legend(loc='upper left',fontsize=legend_meret)

title("Oscillating magnetization $M(B)$ at $T$",fontsize=20)

xticks(fontsize=xyticks_meret)
yticks(fontsize=xyticks_meret)

xlabel(r'$\frac{\hbar \omega_c}{E_\mathrm{F}} \sim B$',fontsize=xylabel_meret)
ylabel(r'$M(B)$',fontsize=xylabel_meret)

xlim(0,B_max*1.)
#\ylim(0,7)

grid();

#savefig('Fig_mu_B.eps');  # Abra kimentese
T, ge, mcs, lmax =  0.01 2 1 10
In [ ]: